Quantum Breathing Mode of Interacting Particles in a 
One-dimensional Harmonic Trap 
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Extending our previous work, we explore the breathing mode — the [uniform] radial expansion 
and contraction of a spatially confined system. We study the breathing mode across the transition 
from the ideal quantum to the classical regime and confirm that it is not independent of the pair 
interaction strength (coupling parameter). We present the results of time-dependent Hartree-Fock 
simulations for 2 to 20 fermions with Coulomb interaction and show how the quantum breathing 
mode depends on the particle number. We validate the accuracy of our results, comparing them to 
exact Configuration Interaction results for up to 8 particles. 
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I. INTRODUCTION 

Harmonically confined few-particle quantum systems 
and especially their time-dependent properties are an 
important subject of experimental and theoretical re- 
search activities. For example, correlated electrons in 
metal clusters^ or quantum dota^Sl and ultracold Bose 
and Fermi gases in traps or optical lattices^ have 
been investigated in recent years. Particularly, Bose- 
Einstein condensation in low dimensions^ and nonide- 
ality (interaction) effecta^"^, including superfluidity and 
crystallization^, lately raised attention. 

Among these properties the behaviour of the breathing 
mode (BM) attracts special interest, as it is easily excited 
experimentahy 9 and turns out 13 to give information on a 
variety of system properties. The BM describes a radial 
expansion and contraction of a finite system and is char- 
acterized by two frequencies in the general quantum case. 
In our previous worlPSl we have shown for a 2-particle 
system that one of those frequencies changes with the 
system dimensionality, the particle spin and the strength 
of the pair interaction. These results were extended^ to 
4 and 6 particles and to different inverse power law in- 
teraction potentials w(r) oc r~ d with d = 1, 2,3. In this 
paper we present new results for larger particle numbers 
and we show how for different coupling strenghts the cor- 
responding breathing frequencies depend on the particle 
number. For that purpose we present the results of time- 
dependent Hartree-Fock simulations for up to 20 particles 
and compare them to the results of exact CI results for 
up to 8 particles. 



II. THEORY 

A. Time-dependent Schrodinger equation 

We br iefly recall the theoretical background of the 
B M 13 ! 14 ! Generally, a system of N interacting particles 



can be described by the Hamiltonian 

N N 



where 
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U + v(Yi 



(1) 



(2) 



is the single-particle Hamiltonian and w(\vi — r\,|) is a 
binary interaction potential. In this case, the external 
single-particle potential v (rj) is chosen to be harmonic, 
i.e. 



Ufa) = -mtfr 2 



(3) 



v serves as a trapping potential with the trapping fre- 
quency fl. In the following, we concentrate on Coulomb- 
interacting particles with equal masses m and equal 
charges e, e.g. electrons. Thus, the interaction poten- 
tial has the form 

•"" r --^ = j^ < 4 » 

with a = e 2 /(47T£ ). Finally, using the notation r = 
(ri, . . . , rjy), the iV-particle time-dependent Schrodinger 
equation (TDSE) reads 
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For convenience, we introduce the scaled quantities f j = 
Ti/lo and t = fit, so that after omitting the tilde symbol, 
the TDSE can be written in the dimensionless form 
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where fo = (h/ (mQj) ' is the well-known oscillator 
length and 



A = 



malo 



(I) 



is the so-called coupling parameter. Due to the rescal- 
ing, there will only be dimensionless quantities through- 
out this work. For example, lengths, times and energies 
are given in units of /o, f2 -1 and M7, respectively. The 
meaning of A can be interpretated as follows. Defining 
the scale of the potential energy as 



E = l -mn 2 ll 

and the mean interaction energy as 
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one finds 
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Hence, A can be understood as the ratio of the interaction 
energy and the confinement energy. The influence of the 
value of A is described later in Sec. Ill CI The actual exci- 
tation of the breathing mode is realized by a fast switch of 
the trapping potential. For a short period of time £ fr the 
single-particle potential v is completely turned off. As a 
consequence, the time-dependent Hamiltonian takes the 
form 

N 

H{t) =Y," h i+ P(*o - *) + °( t - *°ff )] v{n) 
1=1 

N 

+ J>(|r i -r J -|). (11) 
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During the off-time the particles are driven out of their 
initial state. When the potential is restored, the time- 
dependent expectation value of some quantities start to 
oscillate. In particular for the breathing mode, the os- 
cillation of the single-particle potential energy is dom- 
inated by a beating of two harmonic oscillations. The 
corresponding frequencies will be called u r and lor from 
now on. Both a typical time series of the potential en- 
ergy and the excitation process are demonstrated in Fig. 
[TJ In the next subsections we want to point out some 
algebraically accessible properties of the frequencies and 
their relations to the coupling parameter A. 



B. Separation of center-of-mass motion 

As it has already been showrpSl, the system always pos- 
sesses a universal breathing mode whose frequency has 
the value ujr — 2. In the following, we will derive this 
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FIG. 1. Exemplary time series (N — 2, A = 1) of the potential 
energy {E pot ). The waveform shows a superposition of two 
harmonic oscillations with frequencies uj r and ujr. The inset 
demonstrates the excitation of the breathing mode. 



result and show that this value is independent of the cou- 
pling strength, the system dimensionality and the parti- 
cle number. 

The basic idea is the introduction of center-of-mass 
and relative coordinates for the solution of Eq. ([6]) with 
a product ansatz. The center-of-mass coordinate is given 
by 



R 
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! = 1 



and the set of relative coordinates is given by 

x= (xi, x 2 ,.. . ,Xjv-i) , (13) 

with the definition Xj = r^ + i = r^+i — r^. Thus, ^(r,^ 
in Eq. |6]) is replaced by ^(R,x, t). Now the transfor- 
mation is shown for each term in Eq. wm, starting with 



1 N 



■2^dv\ 



2NdK 2 



For the second term we obtain 
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where r 2 k still has to be expressed in relative coordinates 



x. This can be be done as follows: 
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Finally, the third term takes the form 
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As a result the Hamiltonian can be split in two indepen- 
dent contributions: 
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and 
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Hence, the TDSE takes the form 
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The product ansatz 

tt(R,x,i) = 0(R,t)¥>(x,t) 

yields the independent problems 
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and 
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The center-of-mass problem can be transformed to the 
standard oscillator form 
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where the rescaling R = v AR has been used. Now 
consider an initial state, which can be expressed by 



(R,t = 0) = ^c„</> n (R) 



where <f> n is a solution of 
1 d 2 U 
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The associated energy eigenvalues are well-known: E n — 
n + d/2 for all n € {0,1,2,... }, and the time evolution 
of the state in Eq. (pBl is given by 
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The breathing mode manifests itself in the dynamics of 



• iV 



the quantity r 2 = X)i=i r f- Using center-of-mass and 
relative coordinates, this quantity can be expressed ac- 
cording to Eq. (15). For the first term in Eq. (15) a 



breathing frequency can be obtained by determining the 
expectation value 



(R 2 )(t)=iV- 1 / 2 (R 2 )(i) 
The result for this expression is 



with 
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It can be shown with a reduction to the matrix elements 
of a one-dimensional harmonic oscillator that only the 
cases i = j and i = j ± 2 contribute to the last summa- 
tion in Eq. ( 29 ) (the case i = j does not correspond to an 
oscillation) . The only frequency appearing in the oscilla- 
tion is thus given by lur — 2. As the coupling parameter 
A does not appear in the center-of-mass problem, the 
center-of-mass mode with frequency u>r — 2 is present 
for all couplings. In summary, we have shown that the 
frequency uir is universal, but its amplitude tends to van- 
ish for large particle numbers since it is proportional to 

JV-l/2. 



C. Influence of coupling parameter and limiting 



As we have seen in the last subsection, the introduc- 
tion of relative and center-of-mass coordinates has led 
to a separation ansatz. The center-of-mass Hamiltonian 
H-r. yields a breathing frequency ujr = 2. It has already 
been mentioned that the breathing mode also exhibits 
another frequency uj r . The properties of this frequency 
are an interesting subject of numerical analysis. Only in 
two limiting cases the values of u> r are known from alge- 
braic calculations. In the pure quantum limit, A = 0, the 
particles are completely uncorrelated. As the interaction 
term in the Hamiltonian is cancelled out, a degenerate 
frequency u r = ujr = 2 occurs. On the contrary, in the 
classical limit A — > oo the frequency uj r has the value 
V / Jp2"2]. In both cases the frequency does not depend 



on the particle number or the dimensionality of the sys- 
tem. For arbitrary coupling parameters the values of u r 
are expected to be in the interval [v3, 2] . To clarify the 
influnce of the system parameters such as the coupling 
parameter and the particle number on uj t is one of the 
main goals of our investigation. 



III. SIMULATION METHODS 

Whereas the solution of the time-dependent 
Schrodinger equation could only be performed for 
2 particles, the frequencies for up to 8 particles are 
still accessible through exact Configuration Interaction 
calculations. These results are used to support the 
accuracy of time-dependent Hartree-Fock calculations, 
which have been performed for even higher particle 
numbers (up to 20). 

Due to the high computational effort, we restrict our- 
selves to the solution of a 1-dimensional system. Never- 
theless such a system can be regarded as a basic theoret- 
ical model which requires a deepened understanding. In 
order to suppress spin effects, only spin-polarized systems 
are investigated. Before the breathing mode is excited, 
the system is in the energetically lowest anti-symmetric 
state. 

In the following, we give a brief discussion of the em- 
ployed methods and present their numerical results. It 
shall already be mentioned here that in order to avoid 
divergencies in the interaction potential all methods use 



a regularized Coulomb potential A/|r,; — Yj 
k is a small finite cut-off parameter. 



c 2 |, where 



A. Time-dependent Schrodinger equation 



In our previous worH^we determined oj r (X) for N = 2 
in the whole range A = . . . oo. These values are the ba- 
sis for the comparison with other methods. Our TDSE 
results were obtained by solving the time-dependent 
Schrodinger equation with two different methods. On the 
one hand, a standard grid-based Crank-Nicolson scheme 
was used, and on the other hand the wave function was 
expanded into a set of basis functions (oscillator eigen- 
functions) . The results confirm the values of the breath- 
ing frequencies in the limiting cases (A = and A = oo) 
and yield a continuous function w r (A) for all other cou- 
plings in between. 



B. Configuration Interaction 

Configuration Interaction (CI) is another method for 
obtaining numerically exact results. The basic idea of 
CI is to expand the wave function in a complete set of 
Slater determinants, which in turn are constructed with a 
complete single-particle basis. It has already been stated 
in Eq. ( 11 ) that the excitation of the breathing mode is 



realized by a fast Heaviside-type switch-off of the single- 
particle confinement potential. Assuming the excitation 
to be infinitely short, the Hamiltonian is given by 



(31) 



H(t)=H + r)Hi5(t), 



where the time- independent part H is that of Eq. ([lj) . If 
Ho is diagonalized by the eigenfunctions \^f n ) with eigen- 
values E n , the application of perturbation theory yields 
that the expectation value of an arbitrary observable can 
be calculated by 
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where Cij are time-independent coefficients. As a conse- 
quence of this relation, the oscillation of the expectation 
value is restricted to frequencies Wy = \Ei — Ej\. Instead 
of time-propagating the solution of the Schrodinger equa- 
tion, we can use this result and extract the breathing fre- 
quencies from the eigenvalues of Hq with relatively little 
computational effort. However, this method is only appli- 
cable for small particle numbers because of the strongly 
increasing size of the required basis sets for higher par- 
ticle numbers. All results were produced with a basis of 
oscillator functions. Especially for weak couplings this 
basis set is well adjusted to the physical problem and the 
number of basis functions can be kept low. Just like the 
TDSE results the CI results can be used as a benchmark 
for the accuracy of the time-dependent Hartree-Fock re- 
sults. 



C. Time-dependent Hartree-Fock 

For larger particle numbers (N > 7) approximation 
methods have to be employed, because the basis size for 
reasonable CI results dramatically increases. In the fol- 
lowing, the interaction of the system is approximated on 
the mean-field level. The method of choice is the use of 
Nonequilibrium Green's functions G < (1,2), where 1 = 
(ri,ti). Employing the Hartree-Fock (HF) approxima- 
tion, the Green 's fun ctions obey the Keldysh/Kadanoff- 
Baym equation d 18 * 19 ! 



(id tl -h™(l))G<(l,2) 



, (33) 

^°(1) + ^ HF [G<](1), (34) 



where h°(l) is the single-particle part of the Hamiltonian 
(fl]). In such a reduced quantum statistical description 
of the system all one-particle quantities (and some N- 
particle quantities, e.g. the total energy) can be derived 
from the one-particle density matrb<P2 

F 1 (r 1 ,r 2 ,t) = -i lim G < (r 1 , t, r 2 , t') , (35) 
t'->t+ 

where the limit is taken from above. The advantage of the 
quantum statistical description is the fact that the par- 
ticle number is only a parameter and does not determine 



the size of a whole system of equations as in an equiva- 
lent wave-function-based method. The Green's function 
is expanded in FE-DVR (finit e-elem ent discrete variable 
representation) basis function d 21 * 22 * with up to ^400 ba- 
sis functions. The numerical procedure is to calculate 
the ground state (with inverse temperature /3 — > oo), 
switch off and on the confinement and continue propa- 
gating the Green's function in time. In each step the 
expectation values (E pot ) are calculated and saved. Fi- 
nally, one can analyze the spectrum of the time series 
and extract the breathing frequency out of it. Compared 
to methods based on perturbation theory, which only 
involve the ground state energies, this method is quite 
time consuming. Since the computation of just a single 
frequency can last more than one day on one CPU, the 
paramaters have to be chosen carefully in order to guar- 
antee converged results while keeping the durations of 
the computations acceptable. 



IV. RESULTS 

As mentioned before, we want to concentrate on pre- 
senting the results of our time-dependent Hartree-Fock 
(TDHF) calculations and show the dependency of uj t on 
the particle number. In this work we only investigate the 
cases A = 0.1, 0.3 and 1. Larger values of A would re- 
quire to go beyond Hartree-Fock, which demands a very 
high computational effort. Before we show our results, we 
start by explaining some important aspects concerning 
the spectral determination of the breathing frequencies. 
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FIG. 2. Comparison of the spectra for 2, 10 and 18 particles 
at different coupling strengths A. The left peaks represent the 
relative mode with frequency ui r . Peaks with slightly higher 
frequencies than the center-of-mass frequency uje, — 2 can 
be found in case A = 1. For A = 0.1 these peaks and the 
center-of-mass peaks might overlap. 
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FIG. 3. Spectrum corresponding to time series in Fig. [II 
(TDHF) in comparison with the spectrum of the solution of 
the Schrodinger equation (TDSE) for the same configuration. 
The two breathing frequencies ui r and uir can be identified 
in both spectra. However, the TDHF spectrum exhibits nu- 
merous additional peaks. One of those (unphysical) peaks is 
near u) = 2 and has a high amplitude. Aside from that both 
spectra show higher harmonics, which, however, do not agree 
with each other. 



A. Spectral Analysis 

The spectra of all TDHF results were calculated from 
time series which have the length of at least i prop = 2000 
(in units f2 -1 ). Since the resolution in the frequency 
space is limited by the size of t pr0 p, we applied spline in- 
terpolations to the spectra in order to achieve a higher 
accuracy for the frequency values. Besides, each spec- 
trum was calculated with a Blackmail window in order 
to uncover peaks with small spectral weights. In Fig. [2] 
the interesting part of the spectra around the breathing 
frequencies is shown for different couplings and particle 
numbers. The peaks corresponding to w,. can be clearly 
identified. However, the identification of the center-of- 
mass peaks corresponding to u)r = 2 is problematic. 
Only in the case A = 1 it is evidently possible to find 
peaks at frequency ui = 2, although they have a very 
small spectral weight and tend to vanish with increasing 
particle number. Surprisingly, peaks with slightly higher 
frequencies than uj — 2 and a rather strong dependency 
on the particle number are observed. As their spectral 
weight is some orders of magnitude larger than that of 
the center-of-mass peaks, one has to be careful not to 
confuse those peaks with the center-of-mass peaks. This 
is important when regarding the spectra in cases A = 0.1 
and A = 0.3, because for these couplings the small peaks 
at to — 2 vanish. It is conceivable that these peaks and 
those slightly higher than w = 2 tend to overlap. 

The question remains why the spectra show such 
strong additional peaks. In Fig. [3] the spectra of both 
TDSE and TDHF calculations are shown. Both spectra 
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FIG. 4. Part of the spectrum around the breathing frequen- 
cies for different particle numbers at coupling strength A = 1. 
While the relative mode frequency io r slightly changes with 
TV (black line), the center-of-mass mode frequency lor (black 
rectangle) remains constant (apart from numerical errors). 
Note the strongly TV-dependent additional peaks on the right, 
which must not be confused with the lor peaks. 



were produced with the same parameters (A = 1,7V = 
2,k 2 = 0.1). As it can be seen, the TDHF spectrum 
contains several additional peaks. It turns out that the 
appearance of (unphysical) frequencies is a typical effect 
of the TDHF approximation. Consequently, it can be 
a challenging task to distinguish between real physical 
frequencies and artefacts of the approximation method. 
However, in both spectra the peaks of higher harmonics 
occur (Only the first higher harmonics can be seen in the 
figure, but the full data even contain higher harmonics). 
They are caused by the finite duration of the initial exci- 
tation. Unfortunately, the values of the higher harmonics 
are not in a good agreement with each other. Taking into 
account that the spectral weight of the higher harmonics 
is quite small, it is presumable that the accuracy of the 
calculations is not high enough to properly represent the 
higher harmonics. In the following, we want to focus on 
the breathing frequency u r and neglect all higher order 
spectral features. 



B. TV-dependency of the breathing frequency 

The TDHF calculations enabled us to obtain the 
breathing frequencies for up to 20 particles. Apart 
from the problem of clearly identifying the center-of-mass 
mode, we found a typical behaviour of the breathing fre- 
quency ui r , which is qualitatively the same for all cou- 
plings. The change of the spectra with increasing par- 
ticle number is illustrated in Fig. |4J It is obvious that 
the breathing frequency uo r decreases to a minimum (for 
5 particles, see also Fig. pi) and afterwards monotoni- 
cally increases. The figure also shows the center-of-mass 
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FIG. 5. Breathing frequency uj r vs. particle number TV for 
coupling parameters A = 0.1, 0.3 and 1 with cut-off parame- 
ter k 2 = 0.1. For A = 1 the frequencies obtained for different 
k 2 are also plotted. For comparison, the TDSE frequencies 
for two particles and the CI results are plotted as well. For 
k 2 — 0.1, the TV-dependencies are qualitatively the same for 
all A. 



peaks, which, as expected from Eq. (28), tend to vanish 
for high particle numbers. Furthermore, the additional 
peaks near u> — 2 exhibit a rather strong TV-dependency, 
but appear to be converged. 

After analyzing the spectra for the couplings A = 0.1, 
0.3 and 1, it turns out that the TV-dependency of u r is 
qualitatively equal for all A. It is typical that the fre- 
quencies attain a minimum for 5 particles, before they 
start steadily growing. For a complete overview the ex- 
plicit values of the breathing frequencies are summarized 
in the graphs of Fig. [5] For comparison, we also show 
the TDSE values for TV = 2 and the CI values for up to 
8 particles as well as the breathing frequencies obtained 
with other cut-off parameters. As expected, the TDHF 
and the CI values are almost the same (see Fig. [5Va)) for 
small A. With increasing A a constant shift between both 
results arises. The CI results confirm that the breathing 
frequency has a minimum, which, however, occurs for 
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FIG. 6. iV-dependent ground state properties of the inter- 
acting system. The inverse coupling parameter a has a mini- 
mum for 7 particles and afterwards increases with N (a). The 
difference of the addition energy continuously decreases and 
possibly reaches the ideal value 0.5 for large particle numbers 
(b). The pointwise deviation of the densities d has a max- 
imum for 4 particles and tends to vanish for large particle 
numbers (c). 



6 particles instead of 5 particles in the case of TDHF. 
Moreover, the variation of k indicates that the results 
are converged in the regime k 2 = 0.01... 0.1, because 
the results just marginally differ from each other. 

For 20 particles the frequency Lo r has not yet reached 
an asymptotic value. There are, however, several indica- 
tors that for N — > oo the frequency w r converges to the 
value uj — 2 characteristic for the ideal quantum limit. In 
the following, we just consider the ground state. This as- 
sumption is justified by the fact that the breathing mode 
is computed in linear response, which should be derivable 
from the ground state properties. 

Firstly, if the system converges to a pure quantum 
state, the interaction energy must become a negligible 
quantity. Assuming that the Hartree-Fock part of the 
total energy represents the interaction energy, it is pos- 
sible to calculate the ratio of the one-particle energy E 1 
to the interaction energy -Erf: 



a = 



(E 1 ) = (g kin ) + (E pot ) 
(Ehf) (Ehf) 



(36) 



This quantity can be interpretated as the inverse of an 



effective coupling parameter. The values of a for up to 
70 particles are shown in Fig. [6]ja). As expected, a 
increases with the particle number. Interestingly, a has 
a minimum for 7 particles. 

Secondly, one can easily show that for an ideal quan- 
tum system the total energy per particle increases by 0.5 
each time a particle is added to the system, i.e. 
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N 
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N 
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N-l 



= 0.5. 



(37) 



Calculating this addition energy for the interacting sys- 
tem, we get a result which is presented in Fig. [61(b). 
The values slowly converge against the ideal value 0.5, 
although still larger values of N would be required to 
prove this limit. 

A third indicator is the particle density n(r) (or in the 
one-dimensional system n(x)). We want to show that for 
N — >• co the density converges to the density of the ideal 
system. An appropriate quantity for that purpose is the 
pointwise squared deviation of the normalized densities, 
i.e. 



AM 



<ZJ (""ideal (^i) ~ "-interacting (Xj)) ■ (38) 



The factor 1/N prevents this quantity from diverging and 
gives it a relative character. (Recall the normalization 
J n(x)dx = N.) The values for up to 500 particles are 
shown in Fig. [61(c). As expected, d tends to vanish for 
large particle numbers and noticeably has a maximum 
for 4 particles. In order to demonstrate the change in 
the densities, the normalized densities of the ideal quan- 
tum system and the interacting system are plotted for 
2 and 30 particles in Fig. [7] Although all of the indi- 
cators are still not able to prove that for N — > oo the 
system converges to an ideal quantum system, they give 
an evident hint to this hypothesis. An explanation for 
the transition to a pure quantum system might be that 
for large particle numbers the system reaches an exten- 
sion which is dominated by the trap potential (oc rf). 
This is a consequence of the Pauli exclusion principle for 
fermions. Moreover, it is apparent that the quantities a 
and d have an extremum near the minimum of the breath- 
ing frequency. Yet it is not possible to clearly reconstruct 
that mimimum from the ground state. 



C. Conclusions 

We showed that the breathing frequency u> r bears an 
A-dependency which is similar for the coupling strengths 
A = 0.1, 0.3 and 1. Presumably, the same behaviour can 
be observed for arbitrary finite couplings. We found that 
the breathing frequencies decrease to a minimum for 5 
particles (N = 6, for CI) and afterwards monotonically 
increase. 

The origin of the minimum of uj r is the competition 
between two restoring forces: first, the Coulomb repul- 
sion between the particles which is strongest for A — > oo 
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FIG. 7. Comparison of the normalized densities of an ideal 
quantum system and an interacting system (A = 0.1) for 2 and 
30 particles. The deviation from the ideal system vanishes for 
large particle numbers. 



(point charge) and which favours u> r — \/3. Second, the 
kinetic energy of a quantum system that is maximal for 
A = 0. The minimum appears where the ratio of the two, 
1/a has its maximum. It may be expected that the oc- 
curence of this maximum is related to particularly stable 
cluster configurations. 



Although we did not investigate more than 20 parti- 
cles, we found some hints that for N — > oo the breathing 
frequency converges to the frequency uj — 2 of an ideal 
quantum system. In contrast to the increase of tu r with 
N at finite A, for A — > oo, the frequency is independent of 
N and equals v3. Thus, presumably, there exists a crit- 
ical value of A where the crossover between the two be- 
haviours occurs. It remains a subject of further research 
activities to confirm these hypotheses. Furthermore, it 
would be interesting to extend this analysis to higher di- 
mensions, other interactions (e.g. dipole interaction) and 
other spin properties. 
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